SGI Hot Mix 17/Hot Mix 17.iso/HM17_SGI/research/lib/amoeba.pro

Text File | 1997-07-08 | 7KB | 197 lines

; $Id: amoeba.pro,v 1.4 1997/04/11 19:57:21 dave Exp $
; Copyright (c) 1996-1997, Research Systems, Inc.  All rights reserved.
;    Unauthorized reproduction prohibited.


Function amotry, p, y, psum, func, ihi, fac
; Extrapolates by a factor fac through the face of the simplex, across
; from the high point, tries it and replaces the high point if the new
; point is better.

  fac1 = (1.0 - fac) / n_elements(psum)
  fac2 = fac1  - fac
  ptry = psum * fac1 - p[*,ihi] * fac2
  ytry = call_function(func, ptry)  ;Eval fcn at trial point
  if ytry lt y[ihi] then begin    ;If its better than highest, replace highest
    y[ihi] = ytry
    psum = psum + ptry - p[*,ihi]
    p[0,ihi] = ptry
    endif
  return, ytry
end


Function Amoeba, ftol, FUNCTION_NAME=func, FUNCTION_VALUE=y, $
    NCALLS = ncalls, NMAX = nmax, P0 = p0, SCALE=scale, SIMPLEX=p
; The Numerical Recipes procedure Amoeba, with some embellishments.
;
; Reference: Numerical Recipes, 2nd Edition, Page 411.
;  P = fltarr(ndim, ndim+1)
;+
; NAME:
;    AMOEBA
;
; PURPOSE:
;    Multidimensional minimization of a function FUNC(X), where
;    X is an N-dimensional vector, using the downhill simplex
;    method of Nelder and Mead, 1965, Computer Journal, Vol 7, pp 308-313.
;
;    This routine is based on the AMOEBA routine, Numerical
;    Recipes in C: The Art of Scientific Computing (Second Edition), Page
;    411, and is used by permission.
;
; CATEGORY:
;    Function minimization/maximization. Simplex method.
;
; CALLING SEQUENCE:
;    Result = AMOEBA(Ftol, ....)
; INPUTS:
;    FTOL:  the fractional tolerance to be achieved in the function
;    value.  e.g. the fractional decrease in the function value in the
;    terminating step.  This should never be less than the
;    machine's single or double precision.
; KEYWORD PARAMETERS:
;    FUNCTION_NAME: a string containing the name of the function to
;    be minimized.  If omitted, the function FUNC is minimized.
;    This function must accept an Ndim vector as its only parameter and
;    return a scalar single or double precision floating point value as its
;    result.
;    FUNCTION_VALUE: (output) on exit, an Ndim+1 element vector
;    containing the function values at the simplex points.  The first
;    element contains the function minimum.
;    NCALLS: (output) the of times the function was evaluated.
;    NMAX: the maximum number of function evaluations allowed
;    before terminating.  Default = 5000.
;    P0: Initial starting point, an Ndim element vector.  The starting
;    point must be specified using either the keyword SIMPLEX, or P0 and
;    SCALE.  P0 may be either single or double precision floating.
;    For example, in a 3-dimensional problem, if the initial guess
;    is the point [0,0,0], and it is known that the function's
;    minimum value occurs in the interval: -10 <
;    X(0) < 10, -100 < X(1) < 100, -200 < X(2) < 200, specify: P0=[0,0,0],
;    SCALE=[10, 100, 200].
;    SCALE: a scalar or Ndim element vector contaiing the problem's
;    characteristic length scale for each dimension.
;    SCALE is used with P0 to form an initial (Ndim+1) point simplex.
;    If all dimensions have the same    scale, specify a scalar.
;    SIMPLEX: (output and/or optional input) On input, if P0 and SCALE
;    are not set, SIMPLEX contains the Ndim+1 vertices, each of
;    Ndim elements, of starting simplex, in either single or double
;    precision floating point, in an (Ndim, Ndim+1) array. On output,
;    SIMPLEX contains the simplex, of dimensions (Ndim, Ndim+1), enclosing
;    the function minimum.  The first point, Simplex(*,0), corresponds to
;    the function's minimum.
;
; OUTPUTS:
;   Result: If the minimum is found, an Ndim vector, corresponding to
;    the Function's minimum value is returned.  If a function minimum
;    within the given tolerance, is NOT found in the given number of
;    evaluations, a scalar value of -1 is returned.
;
; COMMON BLOCKS:
;    None.
;
; SIDE EFFECTS:
;    None.
;
; PROCEDURE:
;    This procedure implements the Simplex method, described in
;    Numerical Recipes, Section 10.4.  See also the POWELL procedure.
;
;    Advantages:  requires only function evaluations, not
;    derivatives, may be more reliable than the POWELL method.
;    Disadvantages: not as efficient as Powell's method, and usually
;    requires more function evaluations.
;
;    Results are performed in the mode (single or double precision)
;    returned by the user-supplied function.  The mode of the inputs P0,
;    SCALE, or SIMPLEX, should match that returned by the function. The
;    mode of the input vector supplied to the user-written function, is
;    determined by P0, SCALE, or SIMPLEX.
;
; EXAMPLE:
;    Use Amoeba to find the slope and intercept of a straight line fitting
;    a given set of points minimizing the maximum error:
;
;    The function to be minimized returns the maximum error,
;    given p(0) = intercept, and p(1) = slope:
; FUNCTION FUNC, p
; COMMON FUNC_XY, x, y
; RETURN, MAX(ABS(y - (p(0) + p(1) * x)))
; END
;
;    Put the data points into a common block so they are accessible to the
;    function:
; COMMON FUNC_XY, x, y
;    Define the data points:
;   x = findgen(17)*5
;   y = [ 12.0,  24.3,  39.6,  51.0,  66.5,  78.4,  92.7, 107.8, 120.0, $
;        135.5, 147.5, 161.0, 175.4, 187.4, 202.5, 215.4, 229.9]
;
;    Call the function.  Fractional tolerance = 1 part in 10^5,
;    Initial guess = [0,0], and the minimum should be found within
;    a distance of 100 of that point:
;   r = AMOEBA(1.0e-5, SCALE=1.0e2, P0 = [0, 0], FUNCTION_VALUE=fval)
;
;    Check for convergence:
;   if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
;    Print results.
;   print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
;Intercept, Slope:      11.4100      2.72800
;Function value:       1.33000
;
; MODIFICATION HISTORY:
;    DMS, May, 1996.    Written.
;-

if keyword_set(scale) then begin  ;If set, then p0 is initial starting pnt
   ndim = n_elements(p0)
   p = p0 # replicate(1.0, ndim+1)
   for i=0, ndim-1 do p[i,i+1] = p0[i] + scale[i < (n_elements(scale)-1)]
endif

s = size(p)
if s[0] ne 2 then message, 'Either (SCALE,P0) or SIMPLEX must be initialized'
ndim = s[1]            ;Dimensionality of simplex
mpts = ndim+1            ;# of points in simplex
if n_elements(func) eq 0 then func = 'FUNC'
if n_elements(nmax) eq 0 then nmax = 5000L

y = replicate(call_function(func, p[*,0]), mpts)  ;Init Y to proper type
for i=1, ndim do y[i] = call_function(func, p[*,i]) ;Fill in rest of the vals
ncalls = 0L
psum = total(p,2)

while ncalls le nmax do begin        ;Each iteration
  s = sort(y)
  ilo = s[0]        ;Lowest point
  ihi = s[ndim]        ;Highest point
  inhi = s[ndim-1]    ;Next highest point
  rtol = 2.0 * abs(y[ihi]-y[ilo])/(abs(y[ihi])+abs(y[ilo]))

  if rtol lt ftol then begin        ;Done?
    t = y[0] & y[0] = y[ilo] & y[ilo] = t    ;Sort so fcn min is 0th elem
    t = p[*,ilo] & p[*,ilo] = p[*,0] & p[*,0] = t
    return, t        ;params for fcn min
    endif

  ncalls = ncalls + 2
  ytry = amotry(p, y, psum, func, ihi, -1.0)
  if ytry le y[ilo] then ytry = amotry(p,y,psum, func, ihi, 2.0) $
  else if ytry ge y[inhi] then begin
    ysave = y[ihi]
    ytry = amotry(p,y,psum,func, ihi, 0.5)
    if ytry ge ysave then begin
    for i=0, ndim do if i ne ilo then begin
      psum = 0.5 * (p[*,i] + p[*,ilo])
      p[*,i] = psum
      y[i] = call_function(func, psum)
      endif
    ncalls = ncalls + ndim
    psum = total(p,2)
    endif        ;ytry ge ysave
    endif else ncalls = ncalls  - 1
  endwhile
  return, -1        ;Here, the function failed to converge.
end